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ABSTRACT 

Aims. To study turbulence in the Orion Molecular Cloud (OMC1) by comparing observed and simulated characteristics of the gas motions. 
Methods. Using a dataset of vibrationally excited Ha emission in OMC1 containing radial velocity and brightness which covers scales 
from 70 AU to 30000 AU, we present the transversal structure functions and the scaling of the structure functions with their order. These are 
compared with the predictions of two-dimensional projections of simulations of supersonic hydrodynamic turbulence. 

Results. The structure functions of OMC1 are not well represented by power laws, but show clear deviations below 2000 AU. However, 
using the technique of extended self-similarity, power laws are recovered at scales down to 160 AU. The scaling of the higher order structure 
functions with order deviates from the standard scaling for supersonic turbulence. This is explained as a selection effect of preferentially 
observing the shocked part of the gas and the scaling can be reproduced using line-of-sight integrated velocity data from subsets of 
supersonic turbulence simulations. These subsets select regions of strong flow convergence and high density associated with shock structure. 
Deviations of the structure functions in OMC1 from power laws cannot however be reproduced in simulations and remains an outstanding issue. 
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1. Introduction 


Turbulence plays a central role in star-forming molecular 
clouds, acting both to support the clouds globally and to cre¬ 
ate local clumps and density enhancements that can undergo 
gravitational collapse. Simulations have shown that this latter 
process, known as turbulent fragmentation, may directly deter¬ 
mine the initial mass function (IMF). Insight into the effects 
of turbulence on molecular clouds is thus essential for under¬ 
standing the mechanisms of star formation. Such insight can 
only be gained by a close interplay between observations and 
simulations. 

The characterization of turbulence may be 
achieved by statistical methods. Several techniques. 
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> 20021) . have previously been used to characterize the structure 
of brightness or velocity in molecular clouds. Comparisons 
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those in lGustaftsonetal are based on CO observations, 

tracing relatively cool and low density gas. Data are limited in 
spatial resolution and can only be used to address the physics 
at scales larger than roughly 0.03 pc (6000 AU). 


In the present work we use IR K-band observations of vi¬ 
brationally excited H 2 in the Orion Molecular Cloud (OMC1) 
to make a first comparison between observations and hydro- 
dynamical simulations at the scales where individual stars are 
forming. In the region observed the H 2 emission is optically 
thin in the sense that it is not self-absorbed. There will be 
some obscuration by dust dRosenthal et al.112000 ). whose spa- 
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tially differential nature is not known and is ignored here. The 
observations cover scales from 70 AU to 3-10 4 AU. OMC1 is 
the archetypal massive star-forming region and the best stud¬ 
ied region in the sky. OMC1 is highly active with widespread 
on-going star formation, exemplified by the presence of proto¬ 
stars, outflows and larger scale flows (see iNissen et al.il2005t 
and r eferences therein). In a previous paper IGustafsson et alJ 
2006), using the same observational data as in the present work, 
we quantified the nature of turbulence in OMC1 by calculat¬ 
ing size-linewidth relations, probability distribution functions 
and structure functions. It was shown that the turbulence at the 
small scales covered in these data generally follows the trends 
observed in CO data at larger scales. However, analysis also 
showed clear deviations from the fractal scaling observed at 
larger scales. These deviations could be ascribed to the pres¬ 
ence of star formation and associated structures such as circum- 
stellar disks. Here we use the structure functions of the radial 
velocities and the scaling of the structure function exponents 
from our observational data to compare with a numerical sim¬ 
ulation of supersonic, compressible, hydrodynamic turbulence. 
The scaling of structure functions has earlier been analysed in 
IPadoan et al. ( 2003 ) where the column densities of 13 CO were 
used. 

The structure function of order p of the velocity vector u is 
defined here as 

S p (r) = (|[it(:r) — u(x — r)] ■ e\ p ) oc r^ p (1) 

where e is a unit vector parallel (longitudinal structure func¬ 
tion) or perpendicular (transversal structure function) to the 
vector r, and r = |r|. The average is taken over all spatial 
positions x. The modulus sign in our definition 0 is adopted 
to improve the statistics for odd moments. In our case the data 
consist of projected radial velocities. We measure differences 
in radial velocity across the plane of the sky and we are there¬ 
fore dealing with transversal structure functions. The structure 
functions of fully developed turbulent fields are known to fol¬ 
low power laws in the inertial range, 77 <C r <C T, where 77 
is the dissipation scale and L is the integral scale. The set of 
scaling expon ents, ( p in Eq. 0 , can therefore be determined 
dFrischll 1995i) . The scaling exponents are expected to be char¬ 
acteristic of the turbulence involved and universal for all scales 
and Reynolds numbers. The transversal and longitudinal struc¬ 
ture functions are anticipated to have the same scaling in the in¬ 
finite Reynolds number limit. This may not however be the case 
at moderate Reynolds numbers, where it has been found that 
Cp.iong > Cp.trans for p > 3 in inco mpres sible hydrodynamical 
experiments and simulations dlCerr et al .1120011 and references 
therein). 

iKolmogorovI (1941) found from the energy conservation 
law in incompressible, isotropic and homogeneous turbulence 
that ( 3=1 • However, iDubrullel dl994 suggested that ratios of 
scaling exponents, say ( P /C 3 > are inherently universal, while 
the individual scaling exponents_may not be universal them¬ 
selves. In this connection lFrick et al . ( 1995 ) showed in the con¬ 
text of cascade models that one may have £3 ^ 1 and yet re¬ 
cover scaling laws for the structure functions in good agree¬ 
ment with the She-Leveque model (see below) for the ratio 
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IShe & Leveauel J 19941) described the scaling of velocity 
structure functions in incompressible turbulence by: 
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which is c onfirmed by simulat i ons of nearly incompressible 
turbulence (PadoarTet a I J |2004[ Haugen fll, [2004;jl) and by 
experiments ( Anselrnet^taH n^84nRenzietal]n993l) For su¬ 
personic turbulence Boldvre vl d 2002 h obtained, as an extension 
to the She-Leveque model, the scaling: 
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p/3 
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which is confirmed by observations in the Perseus and 
Taurus molecular clouds dPadoan et ak 2003|) and simulations 
dBoldvrev et all2002alhl:l^doan^talj 2004j)/rhis_tvpe of scal¬ 


ing was originally proposed by Politano & Pouauefll995i) for 

magnetohydrodynamic turbulence, where the dissipative struc¬ 
tures are thought to be two-dimensional current sheets. 

In Sect. [2] we describe the observations and calculate the 
structure functions from the observed velocities. We then use 
the method of extended self-similarity (ESS) of iBenzi et al J 
( 1993 ) to find the structure function exponents and show that 
the scaling of these does not represent any known theoretical 
scaling as represented by Eqs. 0 and 0. In Sect. 0 we de¬ 
scribe the simulation. In Sect. l3.ll we calculate the longitudinal 
and transversal structure functions of the 3D simulation and 
show that the scaling of the structure function exponents is sim¬ 
ilar to that of lBoldvrevj (2 002 ) in contrast to the observations. 
In Sect. l3.2l we choose subsets of the simulations which select 
the shocked gas seen in the observations, project the data onto 
2D maps and calculate the structure functions. We show that if 
only strong shocks are included in the subset the scaling of the 
exponents is now similar to the scaling found in the observa¬ 
tions. In Sect. [4] we discuss our results. 


2. Observations 

2.1. Data 

We use the radial velocity map of the BN/KL region of the 
Orion Molecular C l oud (O MC1) of IGustafsson et all J20031 
12004 ; INissen et alJ J2005I) . The data contain both brightness 
and velocity information and were obtained at the CFHT with 
a Fabry-Perot interferometer in conjunction with adaptive op¬ 
tics using the so-called GriF instrument ija enet et a u Eoog) . 
Observations were performed in the NIR K-band by scanning 
the u=l-0 S(l) H 2 emission line at 2. 12177 m. The field of view 
is 36" x 36"and the pixel scale is (y/035 where 1"= 460 AU. 
The dataset consists of four spatial and velocity resolved im¬ 
ages, which are amalgamated into one field of 89"x_67"or 
0.2x0.15 pc for a distance of Orion of 460 pc ( Ballv et alJ 
l200d ). The field is centered approximately on the Becklin- 
Neugebauer (BN) object (05 h 35 m 14H, -05°22'22"9), see 
Fig. [I] The spatial resolution is O'/15 (70 AU). The radial ve¬ 
locity at each spatial position was determined by the peak po¬ 
sition in a lorentzian fit to the velocity profile provided by the 
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Fig. 1 . The velocity field for the full observed region of OMC1. (0,0) is the position of the Becklin-Neugebauer object (BN) 05 h 35 m 14H2, 
—05°22 , 22 , .'9 (J2000). Axes are labelled in arcseconds and colours represent radial velocities in kms . For the white regions no data are 
available. 


Fabry-Perot. Relative velocities are determined with an accu¬ 
racy of between 1 kms -1 (3cr) in the brightest regions and 8- 
9kms -1 in the weakest regions considered here. Systematic 
errors due to mechanical instabilities in the Fabry-Perot may 
occur in establishing velocity differenc es bet ween distant re¬ 
gions. As discussed in lGustafsson et all l2006i ). tests have been 
performed to show that such systematic errors are not present 
in the data to any significant extent. 

The emission of vibrationally excited H 2 observed here 
does not trace the bulk of the gas. Rather it traces hot, dense 
gas, where excitation occurs largely through shock excitation. 
Detailed models Jstorzer & Hollenbach 1999 ) show that the 
maximum brightness in the H 2 v= 1 -0 s(l) line from fluores¬ 
cence in this region of O MC1 does not exceed 1 0-15% of the 
total brightness observed ( iKristensen et a! 120031) . Thus photon 
excitation is a minor contributor to the total brightness. 

2.2. Results from observations 

Since observational data provide only radial velocities in the 
plane of the sky in a 2D projection of the gas motions, we ob¬ 
tain only the transversal structure functions, as described ear¬ 
lier. Furthermore the accuracy of the velocity data in any pixel 
depends on the brightness in that riixel. lGustafsson et al.l ( 2006 ) 
showed that more robust structure functions are obtained when 


the velocity differences are weighted by the brightness. On this 
basis we use a modified definition of the structure functions: 

S p (r) = ( B(x)B(x — r) \v(x) — v{x — r)| p ). (4) 

Here v is the line of sight velocity and the average is extended 
over all spatial positions x and all lags r where r = |r|. B(x) 
is the brightness at position x. We thus weight each velocity 
difference by the product of the brightness of the two spatial 
positions involved, thereby giving more weight to the brightest 
regions which exhibit the highest accuracy in the radial veloc¬ 
ity. 

The third order structure function of OMC1, ,S': 3 (r), is dis¬ 
played in Fig. Eh- It is not well represented by a single power 
law showing a clear deviation around 10 3 AU. This is also ev¬ 
ident from the large variations in the local logarithmic deriva¬ 
tives of S p (r), shown in Fig. 03 for p = 1-5, where the deriva¬ 
tives are evaluated numerically using a three-point formula. 
However. lBenzi et~ak ( 1993 ) discovered that the structure func¬ 
tions can be represented as functions of, say, the third order 
structure function, namely 

S p (r)aS 3 (f) (Cp/C5)ESS . (5) 

This is now known as extended self-similarity (ESS). Even 
if the structure functions of Eq. 0 are not power laws over 
any given range, the functions represented by Eq. 0 nev¬ 
ertheless exhibit good power law behaviour. The scaling in 
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Eq. 0 is generally found to extend over a much larger range 
than for the structure functions of Eq. Q. As emphasized be¬ 
fore, self-similarity, as expressed by Eq. 0 , is believed to be 
more fundamental than the self-similar scaling with respect to 
r feenzi et al 1 ll993l) . In Fig. 0 we have plotted the ratio of 
the logarithmic slopes of S p and S 3 , d\n S p (r)/d\n S 3 (r), for 
p = 1-5. If a range in which good power law scaling is present 
is encountered in the various structure functions, the ratio of 
logarithmic slopes should display plateaus in that range at val¬ 
ues of (Cp/C 3 )ess- FromFig.0 we find that the structure func¬ 
tions for p = 1-5 exhibit a reasonably good scaling range from 
r = 160 AU to r = 7000 AU. This range is marked by the dot¬ 
ted vertical lines in Fig. 0 . The scaling exponents are found 
by fits to Eq. 0 in this range. As an example we show in 
Fig- 01 the extended self-similarity plot of Ss(r) vs. S 3 (r) to¬ 
gether with the best fit yielding the slope (Cp/C 3 )ess = 1-06. 
The dotted lines mark the range of the fit. It is however clear 
from Fig. 0 that the scaling gets poorer when the order p is in¬ 
creased. At p = 5 the plateau is rather poorly defined (see also 
Fig. 01) and therefore we cannot determine a scaling at higher 
orders than p = 5. In Fig. 0 we show the scaling exponents 
(Cp/C 3 )ess vs. p (+ signs) compared to the values predicted by 
the She-Feveque model of incompressible turbulence, Eq. 0 
(dotted line) and the Boldyrev model of supersonic turbulence, 
Eq. 0 (dashed line). The scaling exponents derived from the 
velocity in OMC1 clearly deviate from both the She-Leveque 
and the Boldyrev scaling at p > 4. The OMC1 scaling expo¬ 
nents show signs of becoming constant at (Cp/C 3 )ess ~ 1 or 
even slightly decreasing for p > 4, in contrast to the theoretical 
scalings, which are monotonically increasing. 

This result for velocities of hot, shocked gas in OMC1 at 
scales 70 AU - 3-10 4 AU (3 4-10 ~ 4 pc to 0.15 pc) differs from 
the findings of lPadoan et al . ( i200 3). They found that the density 
fields in the Perseus and Taurus molecular clouds as observed 
in CO follow Boldyrev scaling at scales larger than 0.08 pc. 

Below we will show that the unusual scaling found here 
can be reproduced by numerical simulations of supersonic tur¬ 
bulence when only subsets of the simulations representing the 
shocked regions are considered. 




r / AU 



r / AU 
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3. Simulations 


In order to understand some of the peculiar scalings found 
in the observations we now consider data of supersonic 
isothermal compressible turbulence simulations. Such simu- 
latio ns have been p e rform ed by a number of different groups 
JPassot & Pouaueti Il987t IVazauez-Semadeni eTaTI Il995l 
Padoan et alJ 11998c Klessenl 2000l IV azauez-Semade ni et alJ 


2003( Cho&Lazari anl 2003c Kritsuk&Normanl20041) . Here 

we consider simulations that are most closely related to those 
of Ipfangen et al.l ( l2004h )- except that magnetic fields are 
neglected here. The governing equations are 


8u _ n_, „ 1_ 

— + u ■ Vit = —c: Vlnp + / + - V • r, 

at p 


(6) 


<9 In p 
dt 


+ u ■ V lnp = —V • u, 


(7) 


2.0 

a FS 

^ 1.0 
■Ks' 
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0.0 

0 2 4 6 8 

order p 

Fig. 2. a) Third order structure function of the velocities in OMC1. 
b) Logarithmic derivatives of S p (r) for p = 1-5 c) Ratios of the 
differential slopes of Sp(r) to the slope of the third order structure 
function for p = 1-5. The vertical dotted lines mark the interval in 
which the scaling exponents have been fitted, d) Ss(r) vs. <53 (r). The 
dotted lines mark the range of the fit and the solid line is the best fit 
within that range yielding the logarithmic slope, (Cp/Cs)ess = 1.06. 
e) The ESS scaling exponents (+) OMC1, (dotted line) She-Leveque 
scaling, (dashed line) Boldyrev scaling. 
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Table 1. Parameters of the numerical simulations. Resolution, shock 
viscosity, forcing amplitude, Mach number, run time Ai rvm in terms 
of turnover times Tt urn . 


run 

resolution 

Cshock 

fo 

McLrms 

Atrun/Tturn 

i 

256 s 

2 

2 

3 

270 

2 

512 3 

3 

10 

7-9 

320 

3 

512 3 

2 

10 

8-10 

360 


where T,y = 2puSij + ppStj'V ■ u is the stress tensor 
and S ij = \{uij + Uji) — ^SijUk^k is the rate of strain 
ma trix and commas denote partial differentiation. Following 
Nordlu nd & Galsgaard ( 199 5). we assume // to be proportional 
to the smoothed and broadened positive part of the negative 
divergence of the velocity, i.e. 

p = c s hock /max[(-Vtz)+]\, (8) 

where c s h oc k is the artificial viscosity parameter and the 5 
underneath the max operator indicates that the maximum is 
taken over 5 by 5 by 5 mesh width s (or ’’pixe l s”). Th is is 
also the technique u sed by IPadoan & Nordlundl d2002h and 
lHaugen et alJ J2004bl) . The function / denotes a random forc¬ 
ing function that consists of plane waves, normalized by a di¬ 
mensionless amplitude factor fo that will be varied in the dif¬ 
ferent simulations discussed below (see Appendix A). 

The equations are solved on a periodic mesh of size L 3 , 
where L = 2n/ki is the length of the side of the box and k\ is 
the smallest wave number in the domain. We use the PENCIL 
Code, which is a high-order finite-difference code (sixth order 
in space and third order in time) for solving the compressible 
hydrodynamic equations. 1 

We consider runs with different forcing amplitudes, f 0 , 
leading to different root mean square Mach numbers, Ma rms = 
«rms/Cs where c s is the speed of sound, see Tabled The high 
resolution runs, in rows 2 and 3 of Tabled have been evolved 
for about 40 sound travel times, r soun( j = (c^Aii)” 1 , while the 
low resolution run has been conducted for about 90 T SOU nd- The 
sound travel time can be associated with the turnover time by 
noting that r turn — (^rms ^l) 1 — Tsound /Ma rms* 
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3.1. Results from full 3D simulation 

First, we show that the structure functions of the full 3D sim¬ 
ulation follow the theoretical scaling of Boldvrev (l2002t) . For 
a snapshot of Run 1 at t = 70 T SOU nd (corresponding to t = 
210 rt um ) we have calculated the longitudinal and transversal 
structure functions of the 3D simulation using 

5p,long(0 = “h lj t/, z) U x {X-, U•, 4l) P), (9) 

*5'p,trans(0 — (|tty (x-f* (, t/, 2l) Uy (x, t/, z) P) 

+ (\u z (x + l,y,z) - u z (x,y,z)\ p ), (10) 

as in Eq. 0 . 

In Fig. 0t the third order transversal structure function 
is shown. The logarithmic derivatives of S p .trans(l) (Fig.0>) 

1 http://www.nordita.dk/software/pencil-code 


Fig. 3. a) The third order transversal structure function of Run 1 . b) 
logarithmic derivatives of S p (r) for p = 1-5 (in ascending order), c) 
Ratios of the differential slopes of the transversal structure functions 
of order 1-5 to the slope of the third order structure function, d) The 
ESS scaling exponents of the transversal structure function (+) and 
the longitudinal structure function (diamonds) compared to the She- 
Leveque scaling (dotted line) and the Boldyrev scaling (dashed line). 


show no range of scales where plateaus (good power law scal¬ 
ing) are present. Note, however, that in these simulations £2 
is closer to unity than £ 3 . In Fig. & we have plotted the log¬ 
arithmic Slope Of S p ,trans (0 VS. S^transW for p = 1-5, i.e. 
again using the method of extended self-similarity (ESS). A 
range of good scaling is now seen to be present over most of 
the dynamical range from 10-80 mesh widths. The longitudi- 
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nal structure functions are nearly identical to the transversal 
and are not shown here. The scaling exponents found from fits 
to the structure functions in the interval of 10-80 are plotted in 
Fig. 01 for both the transversal structure functions (+) and the 
longitudinal structure functions (diamonds). Both the transver¬ 
sal and the longitudinal structure functions follow the velocity 
scaling for supersonic turbulence of Boldyrev ( 2002 1. 

3.2. Scaling of subsets of the simulations 

We now study structure functions of subsets of the simulations 
and select subsets that resemble best the physical properties 
of the observational data, that is being composed of preferen- 
tiall y shocked gas. S imilar work has already been carried out by 
iKritsuk & Norman ( 200 41. First, in order to compare the sim¬ 
ulations to the observations, we need to project the simulated 
3D velocity components onto a 2D map of only radial velocity. 
The radial velocity in each spatial position is found by averag¬ 
ing the density weighed ^-component (say) of the velocity over 
the z-range. That is, 

pu z dz j J pdz. (11) 

We have checked that this expression yields the same values of 
velocities as the method adopted in the reduction of the obser¬ 
vational data obtained with GriF. In the observations, the true 
H 2 line profile is convolved with the very much broader in¬ 
strumental Lorentzian profile of the Fabry-Perot interferometer 
and the radial velocity is found from a Lorentzian fit. The same 
procedure has been used on simulated velocity profiles through 
convolution and fitting and it has been found in numerous tests 
that the velocities derived are essentially the same as the cen¬ 
troid velocities obtained via Eq. CD}- In Fig0(top left, marked 
’’all”) the resulting 2D map is shown for Run 1 at t = 50 r sollnc i, 
corresponding to t = 150 Tt urn - 

If the turbulence is homogeneous and isotropic the pro¬ 
jected map should be independent of the projection angle. 
However, since the simulations have limited spatial extent, they 
turn out to show residual anisotropy, in the sense that indepen¬ 
dence of projection angle is not assured. Thus the projection 
map and subsequently the structure functions could depend on 
the projection angle. To minimize such effects we have calcu¬ 
lated the structure functions for a number of random projec¬ 
tion angles. An average of 50 angles was taken for Run 1. The 
higher resolution of Runs 2 and 3 should alleviate the prob¬ 
lem of projection angle and averages over only 3 angles were 
taken in these cases. In the following all structure functions of 
projected maps refer to an average of structure functions. There 
could also be projection effects in the observations, but we have 
no choice but to ignore these. 

The average third order structure function [Eq. 0 with 
p = 3], the logarithmic derivatives of S p (r) forp = 1-5, the 
ratio of the logarithmic slopes of S p (r) and So(r), the extended 
self-similarity plot of S 5 (r) vs. Ss{r), and the velocity scaling 
exponents are also shown in Fig.E] passing down the left col¬ 
umn, for the projected simulation of Run 1. In calculating the 
structure functions of the simulations we use B = 1, that is, no 


2 {x,y) = / 

J Z 


brightness weighting in Eq. 0 since the velocities in the sim¬ 
ulations are free of ‘observational’ errors. It is found that the 
scaling of the structure functions of the projected radial veloc¬ 
ity follows that of Boldyrev, as did the transversal and longitu¬ 
dinal structure functions of the full 3D simulation (Fig.0. 

The problem is to identify the subset of structures in the 
simulations which corresponds to the structures represented in 
our observations. As mentioned above we observe a subset of 
the gas in OMC1 consisting very largely of shocked gas. In 
order to make comparison between observations and simula¬ 
tions it is therefore necessary to extract regions in the simu¬ 
lations where shocks occur. Shocks are generated where fast 
material attempts to overtake slower moving material and ma¬ 
terial shows rapid deceleration, that is, where a negative ve¬ 
locity gradient is present. In simulations, shocked regions can 
thus be distinguished as regions with suitably strong negative 
divergence (V ■ u < 0), that is, convergence. 

Thus, in order to compare with the observations, we choose 
different subsets of the simulations consisting of regions 
with s hocks stronger than a certain degree. Rritsuk & Normanl 
1 2004 1 accomplished this by selecting regions where the den¬ 
sity exceeds a certain threshold value. Here, on the other hand, 
we consider only regions that have stronger convergence than 
a given cut-off value, Dq <0. This is achieved by defining 


D = V ■ u/((V ■ u) 


2\l/2 


( 12 ) 


as the relative velocity divergence, and the selection function s: 

{ 0 £ <-1 

s{x,y,z) = < \ + l£(3-£ 1 2 ) -1 < £ < 1 (13) 

[1 £>1, 

where £ = (Do — D)/w, and w is the width of the selec¬ 
tion function. All points with D > Dq + w, which we seek 
to exclude, have s = 0. The selection function is chosen for its 
smooth variation between 0 and 1. We have used other some¬ 
what different forms of the selection function and found that 
the results do not significantly depend on the specific form. 
Figure El shows an example of a local velocity field where 
shocks are present. The figure shows velocity vectors ( u x , u v ) 
in an xy- plane and contours of shocked regions where s > 0 
for Dq = —6.5 and w = 0.7 in Run 1. 

The shock structure in an xy-slice of Run 3 is shown in 
Fig. E] where regions with large negative values of V • u 
(darker regions) represent shocks. The profiles of u x , u y , and 
u z are shown along a horizontal line on which the presence of 
a shock for example at x = 290 is evident through large veloc¬ 
ity changes in u x , u y , and u z over a range of only ~ 5 mesh 
widths. The higher Mach number of these simulations leads to 
larger velocity differences and stronger shocks compared to the 
simulations of Run l. 2 

The radial velocity in each spatial position (x, y) is now 
found by a modified form of Eq. < 00 : 


u z (x,y) = spu z dz / spdz. 


(14) 


J Z / J Z 

1 Movies of the time evolution of the simulations can be found at: 

http://www.nordita.dk/'brandenb/movies/shockdiss/ 
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Fig. 4. Run 1, t = 50 Tsound- Results for projected maps including all points (1st column) and for subsets with Do = —0.7, —3.6, and —6.5 (as 
marked), using iu = 0.7. Top row: radial velocity maps projected in the ^-direction. Second row: third order structure functions averaged over 
50 projection angles (see text). Third row: logarithmic derivatives of S p (r) for p = 1-5 (in ascending order). Fourth row: ratios of differential 
slopes to £ 3 for orderp = 1-5. The grey shades indicate the ranges over which average values of (£ p /£ 3 )ess are determined. Fifth row: 5s(r) 
vs. 53 (r). The ranges of grey shades correspond to those in the fourth row. Solid lines are best fits within the indicated ranges yielding the 
logarithmic slope, (Cp/C 3 )ess- Bottom row: the radial velocity scaling compared to the She-Leveque scaling (dotted line) and the Boldyrev 
scaling (dashed line). 
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Fig. 5. Vectorplot of ( u x ,u y ) in a section of an xy- plane in Run 1. 
Contours outline shocked regions where s > 0 for the cut-off value 
Do = —6.5, w = 0.7, that is, regions where D < —5.8. 


Returning to Fig. 0 for a snapshot of Run 1 at t = 50 r soun d 
this shows maps projected in the ^-direction as examples of 
subsets with w = 0.7 and Do = —0.7, —3.6, and —6.5. The 
value of D in a region depends on how the velocity changes in 
the vicinity of that region. Large differences in velocity over a 
limited region lead to high values of D. Thus the restrictions 
on D can be associated with typical minimum values of the 
velocity change that must occur in a shocked region for that 
region to be included in the structure function analysis. For ex¬ 
ample the physical interpretation of the restriction D 0 = —3.6 
is that in order for a pixel to be included there must be a ve¬ 
locity gradient in the immediate vicinity of that pixel, such that 
| Am| ~ 3kms _1 over Ar = 10 pixels. An estimate of the 
value of the gradient in physical units can be given by observ¬ 
ing that the size of shocks in the simulations is typically 5 pix¬ 
els. As sum ing a physical shock width of C-shocks of 50 AU 
iLacombe et all2004l) . the value of the required velocity gradi¬ 
ent is ~ 0.03kms^ 1 AU -1 . When Do = —6.5, typical values 
are |Am| ~ 5kms _1 again over A?’ = 10 pixels. These values 
are estimated for Run 1 and are found to be higher in the runs 
with stronger forcing and higher Mach numbers for the same 
value of D. 

The map of Do = —0.7, that is, including all shocked re¬ 
gions, displays sharp filamentary structure, compared with the 
map of all points in the simulation, which has smoother con¬ 
tours (see top row of Fig.0}. Excluding the weakest shocks, 
that is, for Do = —3.6, leads to a more broken up appearance, 
and the filaments are clearly visible. When only the strongest 
shocks are considered. Do = —6.5, the radial velocity map 
consists mostly of sheets and filaments. 

Figure 0] also displays the third order structure functions 
for the three subsets as defined by values of Do, as well as the 
logarithmic derivatives of S p (r), the ratio of the logarithmic 
slopes of S p (r) and 63 (r), the extended self-similarity plot of 
So (r) and So(r), as well as the scaling exponents of the struc¬ 
ture functions of the radial velocity. The structure functions are 
averages of 50 projected maps as described above. The third 
order structure function for the full simulation and the three 
subsets all display good power laws. The leveling off of the 


power laws at lags around 100 mesh widths is an artifact due 
to the finite size of the simulation of size 256. No bumps in the 
structure functions are seen. This is in marked contrast to the 
structure functions of the observations where clear bumps are 
present (see the first panel of Fig. [3. 


The ratios of logarithmic slopes show plateaus over about 
an order of magnitude in scale, especially in the lower or¬ 
der structure functions, p = 1, 2. The scaling exponents 
Qp/Co are found by fits to Eq. (0 in the interval of r where 
din S p {r)/d In S 3 (r) shows the best plateau. The plateau is 
found at r = 5 to 80 mesh widths for the full simulation, at 
r = 10 to 80 mesh widths for Do = —0.7, at r = 2 to 20 
for Dq = —3.6 and at r = 2 to 60 for Do = —6.5. The 
ranges used for fitting are indicated in Fig. 0 with grey shad¬ 
ing. The best power law fits to So(r) vs. Ss(r) in the indicated 
ranges are shown in the fifth row of Fig. 0 The value of the 
slope is indicated. The velocity scaling is close to following 
that of Boldyrev when all points in the simulation are included 
and when only shocked regions with D < 0 ( D 0 = —0.7) are 
included. There is, however, a dramatic change in the scaling 
when the restrictions on the strength of the shocks are made 
tighter. For both Dq = —3.6 and Dq = —6.5 the scaling de¬ 
viates strongly from both She-Leveque and Boldyrev scalings. 
This change in behaviour is associated with a shift in the range 
for which a plateau can be seen. Especially in the last column 
of Fig. 0 the scaling is seen to extend all the way from the 
resolution limit (2 mesh widths) to 60 mesh widths. This can¬ 
not be regarded as regular inertial range scaling because the 
usual dissipative subrange, always present in turbulence sim¬ 
ulations, cannot be distinguished. We associate the apparent 
shift of the ranges with the presence of shock structures which 
become more strongly pronounced as the cutoff value, Dq, is 
moved to more negative values. Thus we are beginning to see 
effects due to the use of artificial viscosity when D becomes 
sufficiently negative. These artificially smoothed shocks may 
resemble C-type shocks that occur in the magnetized interstel- 
lar medium and are known to be a commo n feature in OMC1 
(iGustafsson et al.ll2003t iNissen et al.l2005l) . We also note that 
the nominal dissipation scale without artificial viscosity is just 
1-2 mesh widths, which is still below the artificial dissipation 
scale of the shocks of ~ 5 or more mesh widths. The scaling is 
seen to remain nearly constant for p > 3, resembling the scal¬ 
ing found in the observations of OM C1. This is qualitatively 
similar to the results of lKritsuk & Normanl (2 0041) . who found 
systematically smaller scaling exponents for p > 3 when only 
high density regions are considered. This shows that the un¬ 
usual scaling observed in OMC1 can be the effect of observing 
only the hot, shocked gas, and that hydrodynamical turbulence 
simulations without self-gravity or magnetic fields are able to 
reproduce this. 


Figure |7| shows similar results to Fig. 0] but for a snap¬ 
shot of Run 3 at t = 39 T SOU nd- The subsets corresponding to 
Do = —0.20, —2.0, —7.5, and —9.0 show the same trends as 
seen in Fig.0] The tighter the restrictions on the strength of the 
shocks, the sharper and more broken up the filamentary struc¬ 
tures become. The radial velocity scaling flattens strongly when 
Do = —7.5 and —9.0 (bottom row in Fig.|7}. 
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Fig. 6. Left panel: zy-slice of V ■ u through the box of Run 3 at some arbitrarily chosen value of 2 . The horizontal line indicates the location 
along which the profiles of u x , u y , and u z (in kins -1 ) are shown in the three panels on the right hand side. A strong shock is seen at x = 290. 


The grey-shaded areas in the fourth row of Fig. 0 indicate 
two different regions where scaling can tentatively be noted, 
in contrast to the single regions identified in Fig. 0] As noted 
earlier, the dissipation scale for low velocity gradients is es¬ 
sentially given by 1-2 mesh widths. This is equivalent to the 
case in which no artificial viscosity is introduced; see Eq. (|8}. 
The dissipation scale rises to 5 or more mesh widths in loca¬ 
tions where large velocity gradients are encountered. The two 
different scales identified in Fig.^frows 4 and 5) have dimen¬ 
sions of 2 to >10 mesh widths and 30 to >100 meshwidths. 
The smaller of these ranges covers that associated with artifi¬ 
cial viscosity. The scaling at the lower range will thus be af¬ 
fected by the presence of shocks. This may take place through 
the locally enhanced viscosity embodied in Eq. ©. 

We now consider the scaling associated with the two dif¬ 
ferent regions separately. In Fig. El open triangles refer to the 
larger regions and open squares to the smaller regions. For the 
data in the left column for Dq = —0.20, where all regions 
with a positive value of the convergence are included, that is, 
all shocks, it is seen that both scaling regions (small scales, 2- 
15 mesh widths, and large scales, 20-170 mesh widths) show 
scaling behaviour roughly compatible with Boldyrev scaling. 
When we introduce more restrictive thresholds for Dq different 
behaviour is found. The large scales then begin to follow more 
closely the standard She-Leveque scaling. At the same time, 
the small scales, associated with shocks, show the levelling off 
discussed in connection with Fig. 0 for strongly shocked re¬ 
gions, and as seen in the observations. Fig. Et. We therefore 
are able to explain the unexpected scaling found in the obser¬ 
vations through an inherent selection of shocked regions. 

We note that with the greater energy dissipation associated 
with Run 2, slopes for small and large scales differ as for Run 3 
but less markedly. 


Do = -5.1 _ _ D 0 = -7.5 
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Fig. 8. Run 2, t = 37 r SO und ■ Third order structure functions of pro¬ 
jected maps with w = 0.2 and Do = —5.1 and Do = —7.5. 

3.3. Non-power law behaviour in structure functions 

The structure functions of Run 3 (second row in Fig. 0 are 
all well represented by power laws in contrast to the behaviour 
found in the structure functions of OMC1 (Fig.|2}. In one snap¬ 
shot of Run 2, however, similar bumps to those in the obser¬ 
vations are detected; see Fig. 0 The deviations from power 
laws of the structure functions in Run 2 are found in sub¬ 
sets with high values of the threshold, Dq, in the snapshot at 
t = 37 r soun d. Examples are seen in Fig. [8] for Dq = —5.1 
and Dq = —7.5 where it is clear that the third order struc¬ 
ture functions displays bumps around 20 and 60 mesh widths. 
Other snapshots of Run 2 however do not show this feature. 
The bumps in the structure functions indicate the presence of 
preferred scale sizes in the simulation. This means that even 
if the starting point is a fully isotropic hydrodynamic solution 
of supersonic turbulence, then preferred scales can be encoun¬ 
tered by selecting shocked regions that have a typical filamen¬ 
tary length of some hundred mesh widths (see Fig. |6j. 

4. Conclusion 

We have used observational data of shocked H 2 emission in 
OMC1 to show that structure functions at scales 70 - 3TO 4 AU 
(3.4-10” 4 - 0.15pc) exhibit unusual scaling exponents for 
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Fig. 7. Run 3, t = 39 T SO und. Results for projected maps for subsets with w = 0.20 and Do = —0.20, —2.0, —7.5 and —9.0. Top row: radial 
velocity maps projected in the s-direction. Second row: third order structure functions averaged over 3 projection angles (see text). Third row: 
logarithmic derivatives of S p (r) for p = 1-5 (in ascending order). Fourth row: ratios of differential slopes to £3 for order p = 1-5. The grey 
shades indicate the ranges over which average values of (£ p /£ 3 )ess are determined. Fifth row: Ss(r) vs. Ss(r). The ranges of grey shades 
correspond to those in the fourth row. Solid lines are best fits within the indicated ranges yielding the logarithmic slope, (Cp/Cs)ess. Bottom 
row: the radial velocity scaling compared to the She-Leveque scaling (dotted line) and the Boldyrev scaling (dashed line). A: large scales, □: 
small scales. 
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p > 3. The scaling exponents are nearly constant for p > 3 
and smaller tha n predicted by hoth l.She &. Leveaiiellll 994l) and 
Boldyrev (2002). 


In three simulations we have selected shocked regions by 
imposing requirements on the value of the velocity divergence, 
V • u. In certain important respects the simulations presented 
here are then remarkably successful in reproducing the statis¬ 
tical behaviour observed in OMC1. In other equally important 
areas, they fail. Let us first reiterate the success. 

We have found that by only including shocks that are rel¬ 
atively strong (I)q < —1.5), the unusual scaling exponents of 
the observations are reproduced in the simulations. By contrast, 
a scal ing following that of S he & Leveauel il 1 9941) or lBoldvrevI 
(20021 is found when all points in the simulations are included 
in the data. An explanation for this behaviour is as follows. 
Enhanced energy content at small scales, relative to larger 
scales, implies smaller slopes, that is, smaller values of £ p . Both 
the observational data of OMC1, and some of those subsets of 
the simulations selected only to include shocks, show that the 
values of Q, are reduced for p > 4. Since structure functions of 
high order p are dominated by regions of strong velocity differ¬ 
ences, it follows that the observed excess of small scale energy 
is associated with regions of large velocity differences. These 
are likely to be the regions of strong shocks, as is evidenced by 
the fact that reduced values of Q, are most clearly seen in sub¬ 
sets of the simulations that select the most strongly convergent 
high density regions. 

The present work does not however furnish any explanation 
of why departure from the She-Leveque or Boldyrev scaling 
occurs at the specific value of p > 4. It is possible that the crit¬ 
ical value of p is in some way connected with the physical na¬ 
ture of the shocks, for example the fact that they are smoothed 
in the simulations, mimicking the structure of continuous^ (C- 
) type shocks, rather than jump (J-) type shocks (iFlower et„alj 
feoQ/St and references therein). 

We now turn to the failure of the simulations. The structure 
functions of the observations in OMC1 all deviate from power 
laws and exhibit clear bumps around 10 3 AU, exemplified by 
the third order structure function in Fig. |2] This cannot in gen¬ 
eral be reproduced by the simulations. There appears to be two 
possible explanations for this observed behaviour. The first is 
that the deviation from power laws is due to protostar forma¬ 
tion and associated outflows at a preferred scale. The second is 
that the behaviour is in some way inherent in the nature of the 
turbulence as opposed to the presence of protostars. 

Turning to the first suggestion, the process of star forma¬ 
tion pulls structure of a certain size out of the cascade and 
creates outflows, injecting energy into a turbulent background. 
Gravitational energy and angular momentum is spewed out of 
the star via such outflows and turned into local turbulence, 
hence increasing the overall turbulent content of the gas - 
and restarting the whole cascade process. Such outflows are 
of course not present in the simulations. 

The second suggestion, that the deviation from power law 
is somehow inherent in the nature of the turbulence, requires 
that there is some non-statistical element in this medium which 
is otherwise governed by statistical considerations. This may 
arise through our selection of strong shocks as a subset of the 


whole. We have seen in Fig[8]that traces of bumps in the struc¬ 
ture functions are found in one snapshot of Run 2 when highly 
shocked material is selected. The deviations of the structure 
functions from power laws are not as pronounced in the sim¬ 
ulation as in the observations of OMC1. However this finding 
provides some evidence that part of the explanation for the de¬ 
viations from power laws of the structure functions in OMC1 
arises from the fact that we observe preferentially shocked gas. 
As departures from power law behaviour are only evident in 
a single snapshot and not throughout the simulation at other 
times, this suggestion remains tentative. 

In order to explore the reasons for the departure from 
power law behaviour, more advanced simulations are neces¬ 
sary. These should include self-gravity and energy feedback 
from protostellar zones through outflows and should ultimately 
incorporate ionization and magnetic fields. 
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Appendix A: The forcing function 

For completeness we specify here the forcing function used in 
the present paper. It is defined as jBrandenburgfeOO lh 

f(x,t) = Re{W/ fe(t) exp[ifc(f) • a: + i</>(£)]}, (A.l) 

where x is the position vector. The wavevector k(t) and the 
random phase — n < <p{t) < n change at every time step, so 
f(x. t) is ^-correlated in time. For the time-integrated forcing 
function to be independent of the length of the time step St, the 
normalization factor N has to be proportional to St ^ 2 . On 
dimensional grounds it is chosen to be N = /oC s (|fc|c s /<5f:) 1//2 , 
where /o is a nondimensional forcing amplitude. At each 
timestep we select randomly one of many possible wavevectors 
with length between 1 and 2 times the minimum wavenumber 
in the box, k\. The average wavenumber is k[ = 1.5k’i. We 
force the system with transverse non-helical waves, 

fk = (k x e) /v/fc 2 - (fc • e) 2 , (A.2) 

where e is an arbitrary unit vector not aligned with k; note that 
\fk\ 2 = 1. 
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